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ABSTRACT 

A computer code, FEASBL, is developed to maximize a 
non-linear objective function over a convex feasible region. 
The feasible region is defined by a set of non-linear and 
linear constraints on the variables. FEASBL can solve prob- 
lems involving up to fifty variables with a feasible region 
formed by up to fifty non-linear constraints, and fifty 
linear constraints. FEASBL uses a feasible direction method 
as its solution algorithm. 
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I. INTRODUCTION 



This thesis presents a computer code for solving medium 
sized non-linear programming problems, and describes its use 
at the Naval Postgraduate School computer cent'er. 

The general non-linear programming problem is the optimi 
zation of an objective function over a constrained feasible 
region. By the technique of separable programming, many non 
linear problems can be converted to a linear problem and 
solved using established linear programming methods. Many 
problems, however, cannot be solved by this method. Several 
techniques have been developed for solving non-linear prob- 
lems. Most methods involve following the gradient of the 
objective function to generate a sequence of points that 
converges to the optimal point. The method used in this 
code is Zoutendijk's method of feasible directions [1], The 
method of feasible directions can be considered a large-step 
gradient method as opposed to the small-step gradient and 
search methods. Feasible direction methods can be charac- 
terized by: 

1. Starting at an initial feasible point. 

2. Establishing a feasible direction that will improve 
the objective function. 

3. Proceeding to a new feasible point that has a better 
value of the objective function than the previous point. 

4. Converging to the optimal point, or determining that 
the problem is unbounded. 
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Many procedures can be considered feasible directions methods, 
for example, the simplex method for solving linear program- 
ming problems. 

The method used in this code requires that the user know 
the objective function, constraint functions, all their first 
partial derivatives, and an initial feasible point. If the user 
does not know an initial ' feasible point, this' code can be 
used to solve a preliminary programming problem which will 
provide a feasible point. This problem will be described 
later. This code will solve maximization problems of the 
form : 



Max : 


z = F (x) 




Sub j 


to: G^x) < b i , 


iel r 




T 

a. x < b . , 


ieI L 




d . < x . < c . , 

1 - J ~ 3 


j eJ . 



The functions F (x) , and (x) are non-linear functions which 
are supplied by the user as Fortran FUNCTION routines. The 
sets 1^. and 1^ are the sets of subscripts for the non-linear, 
and linear constraints respectively. If the objective func- 
tion is concave, and the non-linear constraints are convex 
functions, then the problem will be a convex programming 
problem, and the method will converge to the global maximum. 
If the problem is not convex, then the method may converge 
to a local maximum. The code is designed to accommodate 
functions with up to fifty variables. The constraint set 
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can include up to fifty non-linear constraints and fifty 
linear constraints. A complete description of how to use 
this code is provided in Section IV. 
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II. MATHEMATICAL BACKGROUND 



The theoretical foundation for the method of feasible 
directions was developed by Zoutendijk in 1960 [1]. A brief 
discussion of the mathematical background as it applies to 
the method used in this computer code is presented. 

A. DEFINITIONS 

A region R is a convex region if for x^ and X£ > two points 
in R, then the point Xj = ax^ + (1 - a)*^ is in R for all a 
such that 0 < a < 1 . 

Let R be the feasible region, and if x is a member of R 
then a direction s_ is a feasible direction if and only if 
there exists a X' > 0 such that x + Is is in R for all 
Xe(0,X* ) . 

A usable feasible direction is a feasible direction _s 
T 

such that V F(x) s > 0. In other words, a usable feasible 
direction at x is a direction in which the value of the ob- 
jective function increases and along which.it is possible to 
move while remaining in R. 

Let F (x) be a concave function, and let the feasible 
region R be convex. Then if there exists no usable feasible 
direction at x in R, then x must be a maximum of F(x) over 
the region R. To show this assume that x is not a maximum 
and there is no usable feasible direction at x. There exists 
a point £ in R such that F (y) > F(x) and since R is convex 
all the points along the line between and 21 mu st be in R. 
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T 

Further, since F( ) is concave, V x F(x) (^ - x) > 0, thus 
the direction y_ - x is a usable feasible direction. This 
contradiction implies that x must be a maximum of F (x) . 

B . THE METHOD 

The method of feasible directions begins at an initial 
feasible point x° , determines the "best" usable feasible 
direction, moves along that direction to the point x that 
maximizes F (x) along the direction without leaving the fea- 
sible region, repeats this procedure until a point x* is 
found which has no usable direction, or determines that the 
problem is unbounded. The method can thus be viewed as a 
four step procedure. 

1. Determine a usable feasible direction s_ at the 
point x . 

2. If no usable feasible direction exists stop, x is 
the maximum of F (x) over R. 

1c 1 k 1c 

3. Determine a step length A^ such that x = x + Aj,s_ 

k+ 1 k 

is in R and x maximizes F(x) along s_ . 

4. If A^ goes to infinity, then the problem is unbounded, 
otherwise go back to step 1 and repeat. 

The direction finding problem will be discussed first, 
and then the step length problem. 

1 . The Direction Finding Problem 

The direction finding problem is to find the usable 
feasible direction that will result in the greatest apparent 
increase in the objective function. This can be formulated 
as the following linear programming .problem. 



11 



Max: a 

Subj to: a — V F(x) T s < 0 

A 

a + V x G i(x) T s 5 0, iel c (x) 

J s < 0, iel^(x) 

s j > 0, j e J (x) 

Sj < 0, j e J + (x) 

. I s - | < 1, jeJ(x) , 

where I c (x) and I^(x) are the sets of binding non-linear con- 
straints and binding linear constraints, and J (x) and J + (x) 
are the sets of j such that x^ is at its lower and upper 
bound respectively. 

If there exists at least one point that is strictly 
interior to R with respect to the non-linear constraints and 
if the optimal value of a is 0 in the direction finding prob- 
lem, then there is no usable feasible direction and we are 
at the optimum point for the non-linear problem. The last 
constraint in the direction finding problem is a normalization 
constraint. In order to insure that the procedure will con- 
verge the components of s must be bounded, and several nor- 
malizations have been studied by Zoutendijk. These normalizations 
are : 

N 1 : s T s < 1 

N2 : I s . I < 1 

1 r - 
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N3: s. < 1 

3 - 



s. >-l 
3 ~ 



if V F(x) . > 0 
x K - J J 

if V F (x) • < 0 
A J 



N4: o < 1 



N5: V G- (x)^ s + o < b. - G-(x); iel (x) 

x i — — - 1 i — ’ c — 



a- s < b. -a. x; iel. (x) 

““1 “ 1 Li 



T T 




The normalization N1 will generally lead to fewer steps then 
N2, N3, or N4 , but the amount of computation for each step 



and Sj > 0 increases the size of the direction finding prob- 
lem if N1 is used, but not for N2 , N5, and N4. As one goes 
from N2, to N4 the number of steps needed increases, but the 
amount of computation for each step decreases. N5 leads to 
a larger problem, and also does not converge as fast as the 
other normalizations. Normalization N2 was chosen for this 
code as the best compromise between the number of steps and 
the amount of work for each step. 

2 . The Modified Direction Finding Problem 



the maximum, some precautions must be taken to prevent the 
phenomenon of "jamming" in which the procedure in effect 
gets "stuck in a corner" and never finds the real maximum. 
This is accomplished by considering nearby constraints in 
addition to the binding constraints when solving the direc- 
tion finding problem. 




To insure that the procedure actually converges to 
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Define the sets: 



I r (x,e) = {iel r : b. - e < G. (x) < b.} 

1 “X *X 

T 

It (x,e) = (iel, : b. - e < a. x < b.} 

L — ’ J L 1 -—l— -l 

J fx,e) = {jeJ : d. < x. < d- + e} 

' w * J J - J - J 

J + (x,c) = {jeJ : c. - e < x. < c.}. 

- J - J - J 

The direction finding problem is then modified to include the 

constraints in the sets I^(x,e), I^(x,c), J (x,e), and J + (x,e). 

Therefore the direction finding problem becomes: 

Max : a 

Subj to: a - V F(x)^s_ < 0 

+ V x G i^-^ T - 0 » iel c (x >e ) 

T 

_a^ s_ < 0 , iel^(x,c) 

Sj > 0 , jcJ (x,e) 

s. < 0 , jeJ + (x,e) 

|s.| < 1 , j e J (x) . 

Then if the solution to the direction finding problem is such 
that 0 < e, the problem is solved again with e/2. Thus we 
stay away from the "corners" until we get near the optimal 
point . 

3 . The Step Length Problem 

Once the direction s^ has been found, the problem 

is to determine how far to move in that direction. If Ap 

k k 

is defined to be the largest step length A such that x + As_ 
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lc V 

is in R, then the step length will be X^ such that F(x + X^s_ ) 

k k 

is the maximum o£ F(x + Xs_ ) : X < Xp. If X^ = ~, then the 

problem is unbounded, and we can conclude that' the maximum of 

F (x) is infinite. The maximum feasible step Xp is found by 

k k 

a simple search along the direction s^ from the point x . 

Then the optimal step length X^ is determined by finding the 

k k k 

X such that V F(x + s ) s =0. Since F(x) is concave, if 

A ” 

k k T k 

V^Ffx + XpS_ ) _s > 0 , then Xp is used for the step length, 
otherwise X^ can be found by regula falsi [2] . 

4 . Stopping Conditions 

The procedure will terminate whenever the solution 
to the direction finding problem yields a = 0, indicating 
that there are no usable feasible directions. If the solu- 
tion to the step length problem is X^ = °°, then the problem 
is unbounded and the procedure is terminated indicating an 
unbounded solution. In addition, if the step length or the 
percent improvement in the objective function becomes smaller 
than user provided limits, the problem is terminated, and 
the current solution is specified. 



C. THE INITIAL FEASIBLE POINT 

If an initial feasible point is not known, pick an x 1 

such that d. < x! < c., and solve the problem: 

J - J - J 



Max : ~ z \ 

Sub j to: a T x 



Pi z l 



< b. 

i 



id 



L 



d. 

3 



< x . < c . 
‘ J " J 
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where 



z x > 0 



Pi = 



0 ; if aT x* < b. 

1 — - 1 

T x ' -b. ; if aT x ’ >b. 



a 

—l — 



—l — 



An initial feasible point for this linear programming problem 
is x' , = 1. If the feasible region is non-empty, then the 

optimal solution of this problem will be = 0, x = x " . Now 
x" will be feasible with respect to the linear constraints, 
if x" is also feasible with respect to the non-linear con- 
straints, then use x" for the initial feasible point. How- 
ever if x" for the initial feasible point. However if x" is 
hot feasible with respect to the non-linear constraints, 
solve the problem: 



Max : 

Subj to: 



where : 






~ z 2 




G i(x) - 


Pi z 2 5 b i ; ieI 


X\ 

Ctil 


< b^ ; iel 


d . < x. 
J " J 


< c . 
- J 


z 2 i 0 




0 


; if G.(x”) 












G i(x") 



b. ; if G- (x M ) > b. . 
i i — J i 



16 



This non-linear programming problem can be solved by the 
method of feasible directions using x " , = 1 as the initial 

feasible point. If the feasible region of the' original prob- 
lem is non-empty, then the optimal solution to this problem 
will be z 2 = 0, and x. = 21° * N° w the original problem can be 
solved using x° as the initial feasible point. 
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III. THE COMPUTER PROGRAM FEASBL 



This section will describe the computer programming aspects 
of the code, a complete set of user instructions is provided 
in section IV. Several programming languages were considered 
for the program including IBM 360/FORTRAN IV, SIMSCRIPT II. 5, 
and PL/1. FORTRAN IV was chosen because its virtually uni- 
versal availability and widespread familiarity will allow 
easy installation at other computer centers, and facilitate 
user modifications in the future. 

FEASBL is structured as a straight line program, consist- 
ing of several sections; input, direction finding problem 
building and solving, step length determination, and output. 

In deciding the size of the problem that the program would 
be able to handle, the size of the direction finding linear 
programming problem was the controlling factor. To insure 
that round-off errors would not degrade the solution, the 
inverse of the basis matrix is maintained in double precision. 
FEASBL is capable of handling fifty variables, with fifty 
non-linear and fifty linear constraints thus requiring a 
basis matrix that is 152 by 152, which requires over 184K of 
core. Desiring to keep the total size below 300K if pos- 
sible, substantial effort was made to minimize the storage 
requirements throughout the program both by not storing any 
value that could be computed when needed, and using the same 
storage space for multiple purposes. The code now requires 



18 



about 290K of core space in addition to that required for 
user provided functions. 

The information required by the code includes the problem 
description, the initial starting point, and some program 
control variables. The problem description is provided by 
specifying the objective function, non-linear and linear 
constraints, upper and lower bounds on the variables, and the 
requirements vectors. In addition, the first partial deriva- 
tives of the objective function and the non-linear constraint 
functions are required. A prototype code was developed which 
did not require derivatives, but which severely limited the 
form of the functions that could be used. It was felt, how- 
ever, that the greatly increased power and flexibility gained 
by not restricting the form of the functions more than offset 
the requirement to provide the derivatives. These functions 
are provided in the form of Fortran FUNCTION routines, thus 
FEASBL can handle any function the user can express as a 
Fortran routine. The remaining data is input in the Fortran 
NAMELIST format. This method was chosen because the unstruc- 
tured format minimizes keypunch errors. 

Each iteration begins by building the direction finding 
problem. The sets I^Cx^e), I^(x,e), J (x,e) and J + (x,e) are 
determined, and then the direction finding linear programming 
problem is built. Rather then store the entire A matrix, 
only those elements which depend on the gradients of the ob- 
jective function and non-linear constraints in I^,(x,e) referred 
to as the core are stored. The linear programming problem is 
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solved using the revised simplex algorithm [3], and the basis 

inverse is maintained through the product form. At each step 

every non-basic variable is priced, and that variable with 

the most negative z^ - c^ is selected for entry. In the 

event of ties the left-most variable is picked. The column 

vector a., for the entering variable is then generated, and 
— 

used to compute = B • a-^. The departing variable is j 
such that x, /y. , is minimum for all y. , > 0. If there 

DjJ,K J > ^ 

are no negative Zj - Cj's, then the optimal solution has 
been found for the direction finding problem. After an op- 
timal solution has been found, the basis matrix is created 
and inverted directly using the IBM subroutine DMINV, and 
the solution is checked for optimality. If the solution is 
not optimal, then the linear program is continued until an 
optimal solution is reached. When the optimal solution is 
achieved, the direction is available, and the program pro- 
ceeds to find the best step length. This is accomplished by 

k k 

first finding Ap, the largest A such that x + As_ is fea- 

K k 

sible, and then determining A^. < Ap that maximizes F(x + A^s_ ) 
The quantity Ap is found by starting with A = 1 and increas- 
ing A by a factor of 2 until a A^p is found such that 
k k 

x + A. IC s is not feasible. The interval between A vc and 

the largest known feasible step Ap is reduced by checking 

A = MA^rp + Ap) for feasibility until the interval is less 

then the minimum step length (STPMIN) . However, if a feasi- 

k k k 

ble step length A is found such that F (x + As^ ) < F (x ) , 
then it is assumed that the maximum point along the direction 
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lies between these two points. Once Ap has been determined 

A^ is computed by the method of vegula falsi. Since F (x) is 

k k 

concave and is increasing at x along the direction s_ the 

T k 

directional derivative (V F (x) s_) is positive at x and will 

X 

k 

be zero at the maximum along £ . If the directional deriva- 

k k 

tive is negative at x + Ap£ , then the maximum lies between 
the two points. A new A is then selected such that A = 

(bA n - aA ? ) / (b - a) where: a = V F(x^ + A s^) ; b = V F (x^ 

+ ^2— )» an< ^ anc ^ ^2 ^ rac ^ et t ^ ie optimal stop length. If 
the directional derivative is zero at the new point the maxi- 
mum has been found, otherwise we have a new interval with 
positive and negative directional derivatives at the end- 
points, and another A is determined. This continues until 
either the maximum is found, or the interval becomes less 

then the minimum specified. If the directional derivative 

k k k k k 

is positive at x + Ap^ and F (x + Ap£ ) > F (x ), then the 

k 

maximum along s_ is outside the feasible region, and Aj, is 

k k k 

set equal to Ap. However, if F(x + Ap£ ) < F (x ) , then the 

function is not concave, so a step A^ between 0 and Ap is 

found such that the directional derivative is negative, and 

k k k 

it is assumed that F (x) is concave between x and x + AyS_ , 
and the maximum is sought as before. 

Once a usable feasible direction, and step length have 
been determined set x = x + A^s_ and repeat the process 
starting at x^ + k This process is continued until either 
no usable feasible direction can be found or the problem 
becomes unbounded. If the step length or the percent 
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improvement in the objective function become less then user 
supplied limits the program will terminate. The program will 
also stop after a user specified number of iterations. 
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IV. USING FEASBL 



This section will describe how to use FEASBL, including 
required input, possible errors, and the normal output. An 
example is given illustrating different features of the code 

A. INPUT 

The objective function, non-linear constraint functions, 
and their first partial derivatives are input as Fortran 
FUNCTION routines. The linear constraint matrix, require- 
ments vectors, upper and lower bounds on x^ , the initial fea 
sible point, and any other input is read in as data at 
execution . 

1. The FUNCTION Routines 

Four Fortran FUNCTION routines are required, these 
routines provide the objective function, its first partial 
derivatives, the non-linear constraints, and their first 
partial derivatives. Even if the problem has no non-linear 
constraints, all four FUNCTION routines must be provided. 

The four routines are: 



FUNCTION 



Definition 



FX( X ) 

DELFX( NX, X ) 



The objective function. 

The partial derivatives of 
the objective function. 



GIX ( NGI , X ) 



The non-linear constraint 
functions . 



DELGIX (NX ,NGI , X) The partial derivatives of 

the non-linear constraints 
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The arguments of the FUNCTIONS statements are: 



Argument 



Definition 



X 



The vector x. 



NX 

NGI 



Subscript of the partial 
derivative desired, i.e. 
3F (x)/3x NX . 

The non-linear constraint 
number . 



Each routine must include a statement dimensioning X(e.g.; 
REAL*4 X(8), or DIMENSION X(15) ), at least one assignment 
statement with the name of the FUNCTION on the left of the 
equal sign (e.g.; FX = X(l), or DELGIX = 4.0 * X(2)**3 ), 
and at least one RETURN statement. No variation in the names 
of the FUNCTIONS, nor in the order of the arguments is al- 
lowed. The following example illustrates the use of the 
FUNCTION routines. 



Max : 

Subj to: 



F(x) 


= 6x^ 


* 8x 2 


- 


4x? 


+ A 


< 


16 


1 


2 






(4 - 


xp 2 - 


4 <- 


16 


3x. 


+ 5x_ 


< 


15 


1 


2 






2x. 


+ 12x_ 


< 


30 


1 


2 






0 < 


x. < 1 . 


85 






1 - 






2.0 


< x 2 . 








/ 



x 



2 

2 



non-linear constraints 



linear constraints 
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For this problem the following FUNCTION routines are required: 

a. The objective function 

FUNCTION FX( X ) 

REAL *4 X (2 ) 

FX = 6.0 * X(l) + 8.0 * X (2) - X (1) ** 2 - X(2)**2 

RETURN 

END 

b. The partial derivatives of the objective function 

FUNCTION DELFX ( NX, X ) 

REAL * 4 X(2) 

GO TO ( 101, 102 ) , NX 

101 DELFX = 6.0 - 2.0 * X(l) 

RETURN 

102 DELFX = 8.0 - 2.0 * X (2) 

RETURN 

END 

c. The non-linear constraint functions 

FUNCTION GIX ( NGI , X ) 

REAL * 4 X (2) 

GO TO ( 100 , 200 ) , NGI 
100 GIX = 4.0 * X (1 ) **2 + X (2) ** 2 
RETURN 

200 GIX = ( 4.0 - X(l) ) * * 2 + X (2) **2 
RETURN 
END 

d. The partial derivatives of the non-linear con- 
straints 

FUNCTION DELGIX ( NX, NGI, X ) 

REAL* 4 X (2) 

GO TO ( 100 , 200 ) , NGI 

100 GO TO ( 101, 102 ) , NX 

101 DELGIX = 8.0 * X(l) 

RETURN 

102 DELGIX = 2.0 * X(2) 
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RETURN 



200 GO TO ( 201, 202 ) , NX 

201 DELGIX = 2.0 * X(l) - 8.0 
RETURN 

202 DELGIX = 2.0 * X(2) 

RETURN 

END 

The significant feature to note is that when a FUNCTION is 
called by FEASBL specifying a point X, and when necessary 
the partial derivative desired NX, and the constraint number 
NGI, the appropriate value is returned. 

2 . Data Input 

The remainder of the data needed is read in from 
data cards at execution in NAMELIST format. This format was 
chosen to minimize keypunch errors, and provide flexibility 
for the user. The following variables are required: 

Name Description 

IXMAX Number of variables. 

IGCNT Number of non-linear constraints. 

ILCNT Number of linear constraints. 

The following variables are optional, and if not provided by 
the user will take on the default values listed. 



Name 


Default 


Description 


GL ( I , J) 


o 

o 


Linear constraint matrix element 


XMIN(I) 


o 

o 


Lower bound on x^. 


XMAX ( I ) 


= 10 7S 


Upper bound on x^ 


BC ( I ) 


0.0 


The b. for the non-linear con- 
. .r 

s traa nts . 
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Name 


Default 


Description 




blci) 


O 

O 


The 


b. for 
i 


the linear constraints 


X(l) 


o 

o 


The 


initial 


feasible point. 


EPSLON 


0.1 


The 


initial 


anti-jamming constant 


STPMIN 


10' 3 


The 


minimum 


step length. 


PERCNT 


IQ’ 2 


The 


minimum 


objective function 






improvement 


• 


ITMAX 


25 


The 


maximum 


number of iterations. 



The first data card begins with an ampersand (S) in column two 
followed by the word INPUT with no space ( SINPUT' ] , The 
values are then listed as "variable name = FORTRAN constant" 
i.e. IXMAX = 2, GL(1,1) = 3.0, the last data item is followed 
by an ampersand and the word "END" ( SEND ). If more then 
one card is required, each card must begin in column 2 or 
more, and the SEND will appear on the last card only. For 
the above example the data cards would look like: 

SINPUT IXMAX = 2, IGCNT= 2 , ILCNT=2, GL(1,1) = 3.0, GL(1,2) = 5.0 
GL (2 , 1) = 2 . 0 , GL (2 , 2 ) = 12 . 0 , XMAX (1) = 1 . 85 , XMIN(2) = 2.0, 

BC ( 1) = 16 . 0 , BC (2) =16 . 0 , BL(1)=15.0, BL(2)=30.0, ITMAX=10, 
EPSLON=0 .01, STPMIN= 1 . 0E - 04 , X(l) = 1.0, X(2) = 2.0, SEND 
Here the order is unimportant, and all the variables need not 
be listed. In this example the default values for XMIN(l), 
XMAX(2) , and PERCNT will be used since they are not listed 
on the data cards . 

3 . Data Limits 

Input errors will result if the variable inputs are 
outside the following limits: 
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1 < IXMAX < 50 



0 < IGCNT < 50 
0 < ILCNT < 50 

XMIN(I) < XMAX(I) 

To prevent stopping prematurely, EPSLON should be significantly 
larger than STPMIN. 

After all the data is input, the initial point is 
checked for feasibility. If the point is not feasible the 
program is terminated. 

4 . The Program CHECK 

FEASBL requires 290K of core space in addition to 
that required for the FUNCTION routines. Most problems will 
require less then 320K and less then 4 minutes to run. As 
an aid for checking the input data and FUNCTION routines , the 
program CHECK can be used. CHECK requires exactly the same 
input as FEASBL, and performs the same input error checks. 

In addition CHECK will determine the feasibility of the ini- 
tial point and print out the value of the objective function, 
all the constraints, and the value of all the partial deriva- 
tives of the objective function and non-linear constraints. 

By using CHECK the user can quickly verify the input, since 
CHECK can be run normally as a class B or C job, before 
using FEASBL, which is a class J job. 

5 . System Cards 

The JCL cards required to use FEASBL or CHECK depend 
on the way that the FUNCTION routines are input. The routines 
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can be input as either Fortran source decks and compiled as 
part of the job, or the routines can be precompiled and in- 
put as object decks. An example of each methad follows. 

a . FEASBL with Source Decks 

// standard OS JOB card ,TIME= (0 3 , 59) 

// EXEC FORTCLG, REGION. GO=300K 
//FORT . SYSIN DD * 

(FUNCTION source decks) 

/* 

//LINK.NLPLIB DD DSNAME=F0600 . NLPLIB ,UNIT=2 32 1 , 

// VOL=SER=CELOO 1 , DISP=SHR, LABEL= ( , , , IN) 

//LINK. SYSIN DD * 

INCLUDE NLPLIB (FEASBL, DARRAY) 

/* 

//GO. SYSIN DD * 

5 INPUT IXMAX=2 , IGCNT = 2 , ... , $END 
/* 
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b . FEASBL with Object Decks 



// standard OS JOB card ,TIME= (03 , 59) 

// EXEC LGO, REGION. GO=300K 

/ /LINK.NLPLIB DD DSNAME=F060 0 . NLPLIB ,UNIT= 2 32 1 , 
// VOL=SER=CELOO 1 , DISP -SHR , LABEL= ( , , ,IN) 

//LINK. SYSIN DD * 

INCLUDE NLPLIB (FEASBL, DARRAY) 

(FUNCTION object decks) 

/* 

//GO. SYSIN DD * 

§ INPUT IXMAX=2 , IGCNT=2 , ... , X(2) = 2.0, §END 
/* 
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c. CHECK with Source Decks 



// standard OS JOB card 
// EXEC FORTCLG, REGION. GO=100K 
//FORT . SYSIN DD * 

(FUNCTION source decks) 

/* 

//LINK.NLPLIB DD DSNAME = F0600 . NLPLIB ,UNIT = 2 32 1 , 
// VOL=SER=CEL001 ,DISP=SHR, LABEL = (, , , IN) 
//LINK. SYSIN DD * 

INCLUDE NLPLIB (CHECK) 

/* 

//GO. SYSIN DD * 

$INPUT IXMAX=2 , I GCNT= 2 , ... ,X(2) = 2.0, §END 
/* 
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d . CHECK with Object Decks 



// standard OS JOB card 

// EXEC LGO, REGION. GO=100K 

/ /LINK . NLPLIB DD DSNAME=F0600 . NLPLIB, UNIT=2321 , 

// VOL=SER=CEL001,DISP=SHR,LABEL=(, , ,IN) 

//LINK . SYSIN DD * 

INCLUDE NLPLIB (CHECK) 

(FUNCTION objective decks) 

/* 

//GO. SYSIN DD * 

§ INPUT IXMAX=2 , IGCNT= 2 , ... , X(2)=2.0, §END 
/* 

A set of object decks can be produced either after 
checking the FUNCTION routines by executing a FORTCD job, or 
in conjunction with the program CHECK by replacing // EXEC 
FORTCLG, REGION. GO=100K with // EXEC FORTCLGP , PARM . FORT= 
'DECK' , REGION. GO=100K. 

B. SAMPLE PROGRAM 

The FUNCTION routines, input cards, and system cards 
required to use FEASBL with source decks are shown for the 
problem: 

2 2 

Max: F (x) = 6x^ + 8X2 - x^ - X2 

Subj to: 4x^ + x 2 f 16 

(4 - x-^) 2 + x 2 - 
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< 15 



3x^ + 5x^ 

2x^ + 12x2 5 30 

0 < x 1 < 1.85 
2 . 0 < x^ < 00 

A sample of the output of both CHECK and FEASBL for this 
problem is provided. The next two pages are the cards re- 
quired for executing FEASBL with the sample problem. This 
is followed by the output of CHECK and then FEASBL for the 
sample problem. 



33 



noon ooooo ooo non ooo ooo ooooo non ooo noon ooooo ooo ooo 



//SAMPLE01 JOB (1115. 0601, R0Y2) ,' SAMPLE JOB • , T I ME = ( 03 , 5 9 ) 
// EXEC FORTCLG, REGION. GO=300K 
//FORT.SYSIN DD * 

FUNCTION FX( X ) 

THIS FUNCTION PROVIDES THE OBJECTIVE FUNCTION. 
REAL*4 X ( 2 ) 

THIS IS THE OBJECTIVE FUNCTION 

FX = 6.0 * X ( 1 1 + 8.0 * X ( 2 ) - X(l)**2 - X(2)**2 
RETURN 
EN0 



jjcajc I*###* $ # # # # # # jjc 5}: 5}: ajt £ # 

FUNCTION DEL F X( NX, X ) 

THIS FUNCTION PROVIDES THE FIRST PARTIAL DERIVATIVES 
OF THE OBJECTIVE FUNCTION. 

REAL*4 X ( 2 ) 

GO TO ( 101, 102 ) , NX 

THIS IS THE PARTIAL WITH RESPECT TO X(l). 

101 DELFX = 6.0 - 2.0 * X < 1 ) 

RETURN 

THIS IS THE PARTIAL WITH RESPECT TO X ( 2 ) 

102 DELFX = 8.0 - 2.0 * X ( 2 ) 

RETURN 

END 



FUNCTION G I X ( NGI, X ) 

THIS FUNCTION PROVIDES THE NON-LINEAR CONSTRAINTS. 
REAL*4 X ( 2 ) 

GO TO THE PROPER CONSTRAINT. 

GO TO ( 100, 200 ) , NGI 

THIS IS NON-LINEAR CONSTRAINT NUMBER 1. 

100 G I X = 4.0 * X(l)**2 + X(2)**2 
RETURN 

THIS IS NON-LINEAR CONSTRAINT NUMBER 2. 

200 G I X = ( 4.0 - X(l) )**2 + X ( 2 ) **2 
RETURN 
END 






FUNCTION DELG IX ( NX, NGI, X ) 

THIS FUNCTION PROVIDES THE FIRST PARTIAL DERIVATIVES 
OF THE NCN-L I NEAR CONSTRAINT FUNCTIONS. 
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non ooo ooo ooo ooo ooo ooo 



REAL** X ( 2 ) 

GO TO THE CONSTRAINT SPECIFIED. 

GO TO ( 100, 200 ) ,NGI 

THIS IS CONSTRAINT NUMBER ONE, WHICH PARTIAL IS WANTED 

100 GO TO { 101, 102 ) , NX 

THIS IS THE PARTIAL WITH RESPECT TO X(l). 

101 DELGIX = 8.0 * X( 1 ) 

RETURN 

THIS IS THE PARTIAL WITH RESPECT TO X ( 2 J 

102 DELGIX = 2.0 * X(2) 

RETURN 

THIS IS CONSTRAINT NUMBER TWO, WHICH PARTIAL IS WANTED 

200 GO TO ( 201, 202 ), NX 

THIS IS THE PARTIAL WITH RESPECT TO X ( 1 ) 

201 DELGIX = 2.0 * X ( 1 ) - 8.0 
RETURN 

THIS IS THE PARTIAL WITH RESPECT TO X ( 2 ) 

202 DELGIX = 2.0 * X ( 2 ) 

RETURN 

END 



//L INK.NLPL I B DD DSN=F0600. NLPL I B , UN I T = 232 1 , 

// VCL=SER=CEL001 ,DISP=SHR, LA6EL= ( , , , IN) 

//L INK . SYS I N DD * 

INCLUDE NLPL IB( FEAS BL , D ARRAY ) 

//GO.SYSIN DD * 

6INPUT I XM AX =2 , IGCNT = 2, ILCNT=2, GL(1,1)=3.0, GL(1,2)=5.0, 
GL (2,1) =2.0, GL(2, 21 = 12. C, XMAX ( 1 ) - 1 . 8 5 , XMIN(2)-2.0, 

BC ( 1 ) = 1 6 .0 , BC ( 2 ) = 1 6 . 0 , BL(1) = 15.0, BL(2)=30.0, ITMAX=10, 

E PSLON =0.01, STPMIN = 1 .0E-04, X ( 1 ) = 1 . 0 , X(2) = 2.0, SEND 
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SAMPLE OUTPUT FROM CHECK 



NON-LINEAR PROGRAMMING PROBLEM 



INITIAL SET UP: 

NUMBER OF VARIABLES: 2 



NUMBER OF 


NON-LINEAR 


CONSTRAINTS: 


: 2 




NUMBER OF 


LINEAR CONSTRAINTS: 


2 




THE LINEAR 


constraint 


MATRIX IS: 






( l, l); 


3. 00000 


( 1, 


2) ; 


5.00000 


( 2, l); 


2.00000 


( 2, 


2); 


12.0000 



THE BOUNOS ON THE VARIABLES ARE: 



0.0 < X( 1) < I.B5000 

NON-LINEAR CONSTRAINTS REQUIREMENTS VECTOR, BC(I): 
< 1); 16.0000 t 2); 16.0000 

LINEAR CONSTRAINTS REQUIREMENTS VECTOR, BL(I): 

( I); 15.0000 ( 215 30.0000 

THE INITIAL FEASIBLE STARTING POINT IS; 



2.00000 



X( 1); 1.00000 

ANTI JAMMING CONSTANT : 
MINIMUM STEP LENGTH: 
MINIMUM IMPROVEMENT: 
MAXIMUM 4 ITERATIONS: 



,X( 2); 

0. 100000E-01 
0. 1 OOOOOE- 03 
0. 100000E-01 
10 



2.00000 



,X( 



INPUT COMPLETE. INITIAL VALUE OF OBJECTIVE FUNCTION: 



< X( 2) < 



0 • 7237 01 E 7 o 



1.700000E 01 



VALUE OF NON-LINEAR CONSTRAINT NUMBER, 
VALUE OF NON-LINEAR CONSTRAINT NUMBER, 
THE VALUE OF LINEAR CONSTRAINT NUMBER 
THE VALUE OF LINEAR CONSTRAINT NUMBER 



NEXT THE VALUES OF THE PARTIAL DERIVATIVES OF THE OBJECTIVE FUNCTION ANO THE CONSTRAINTS ARE LISTEO. 



1. 


IS: 


B . OOOOOOOE 


00 


2, 


IS : 


1. 3000000E 


OL 


I , 


IS: 


1. 300000E 


01 


2, 


IS: 


2.600000E 


01 


INEAR 


constraints 


HAV 



THE 


VALUE 


OF 


0 EL FX ( 


1, XJ IS: 


A . 


OOOOOOE 00 




THE 


VALUE 


OF 


OELFX( 


2, X) IS: 


4. 


OOOOOOE 00 




THE 


VALUE 


OF 


0 EL G I X ( 


I* 1. XJ 


IS: 


B. OOOOOOE 


00 


THE 


VALUE 


OF 


0 EL G I X ( 


2, 1, X) 


IS: 


4. OOOOOOE 


00 


the 


VALUE 


OF 


OELGIXt 


It 2, X) 


IS: 


-6. OOOOOOE 


00 


THE 


VALUE 


OF 


OELGIX( 


2, 2, X) 


IS: 


4. OOOOOOE 


00 
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SAMPLE OUTPUT FROM FEASBL 



NON-LINEAR RROGRAMMING PROBLEM 



INITIAL SET UP: 

NUMBER OF VARIABLES: 

NUMBER OF NGN-LINEAR CONSTRAINTS: 
NUMBER OF LINEAR CONSTRAINTS: 

THE LINEAR CCNSTRAINT MATRIX IS: 

( 1, l); 3.00000 ( 1, 2) 



2 

2 

2 



5.00000 



( 2 , 1 ); 2.00000 



( 2 * 2 > ; 12.0000 



the bounds cn the variables are: 



0.0 < X( 1) < 1.B5000 ; 

NON-LINEAR CONSTRAINTS REQUIREMENTS VECTOR, BC(I): 

( l); 16.0000 ( 2); I6.0000 

LINEAR CONSTRAINTS REQUIREMENTS VECTOR, BL ( I ) : 

.{ 1); 15.0000 ( 2); 30.0000 

THE INITIAL FEASIBLE STARTING ROINT IS: 



2.00000 



< X( 2) < 



xc l); I.00000 ,x( 2 ); 

ANTI JAMMING CONSTANT : 0.100000E-01 

MINIMUM STER LENGTH: 0.100000E-03 

MINIMUM IMPROVEMENT: 0.100000E-01 

MAXIMUM * ITERATIONS: 10 

INPUT COMPLETE. I NI T I AL VALUE OF OBJECTIVE FUNCTION: 



2.00000 



,X( 



1 . 700000 E 01 



PROBLEM SOLUTION 



1 EK F (X^ N = NUKl " 1? B67500E 01, IMPROVEMENT THIS ITERATION: 11.0294 RERCENT . 



THE RCINT X( ) IS: 
X( I), 1.250000 



X( 2) , 



2.2 50000 



X( 



ITERATION NUMBER: 2 

F ( X ) - I.S22217E 01*. IMPROVEMENT THIS ITERATION: 



THE POINT X( ) IS: 

X( 1), 1.666625 ; XC 2), 2.000024 

CRTIMAL SOLUTION: 

THE CRTIMAL VALUE OF THE OBJECT IVE FUNCTION IS: 
THE PGINT XC ) IS : 

XC 1), 1.666625 ; XC 2), 2.000024 



1.8393 PERCENT. 

; xc 

1.922217E 01 

; xc 



0.723701E 76 
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V. EVALUATION 



To validate the code, and evaluate its performance, 
several non-linear problems were solved with FEASBL. These 
problems were selected to represent the range of problems 
FEASBL was designed to solve. The problems included uncon- 
strained optimization of a concave objective function, con- 
strained optimization of concave objective functions some 
with the unconstrained maxima interior and some exterior 
to the feasible regions. An unbounded problem was also 
solved . 

In each case the optimum solution was obtained in less 
then five seconds of cpu time. These problems were all small 
in order to facilitate manual verification of the solutions. 
Larger problems will take correspondingly more time, but the 
code appears to be quite efficient. In most problems the 
optimal solution was obtained after two iterations. In one 
case, the solution took seven iterations. 

It was found that the value of STPMIN is critical to the 

performance of FEASBL. If STPMIN is set too small, when the 

interval bet\'/een and X Mr . is calculated, the round-off 

F NF 

error may exceed the value of STPMIN. This results in the 
code looping indefinitely trying to reduce the interval beyond 
machine accuracy leading to a program completion due to ex- 
cess time. If this should occur the problem should be re- 
started at the last point with a larger value for STPMIN. 

With step lengths on the order of 1.0 the interval cannot be 
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determined to an accuracy greater then 10 With larger 
step lengths this value will increase. An initial value 

- 2 - 7 

of 10 or 10 is recommended for STPMIN. This value 
can be reduced in subsequent runs of FEASBL if a more pre- 
cise solution is desired. 

The initial value of EPSLON is not as critical as STPMIN. 
EPSLON should be greater than STPMIN but generally not greater 
than 10 . Experience with FEASBL indicates that values on 

the order of 10 or 10 perform well. If a large value is 
used, on the order of 1.0 or greater, and an optimal point 
is found, this should be verified using a smaller value for 
EPSLON, as the code may stop early with large values of 
EPSLON. 

Performance was uneffected by variations in PERCNT. Large 
values may cause the program to terminate sooner. Values less 
then 10 ^ are equivalent to zero because of inherent machine 
accuracy . 

An attempt was made to solve a problem in which all the 
constraints were equality constraints, thus the feasible 
region had no strictly interior points. The program would 
terminate immediately without any significant movement. This 
phenomenon was anticipated since the underlying theory assumes 
that the feasible region contains at least one strictly in- 
terior point to the set of non-linear constraints [1]. 
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